function [GT, CD] = Cragg_Donald(Y, X, Z, ZZ, n, k1, k2)

%Define matrices
ZZr=(eye(n)-X*inv(X'*X)*X')*ZZ;
P=ZZr*inv(ZZr'*ZZr)*ZZr';
Yr=(eye(n)-X*inv(X'*X)*X')*Y;
Mz=(eye(n)-Z*inv(Z'*Z)*Z');
S=(Y'*Mz*Y)/(n-k1-k2);

%Compute Cragg-Donald statistic 
GT=(S^(-0.5))'*Yr'*P*Yr*(S^(-0.5))/k2;
CD=min(eig(GT));